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ABSTRACT 


Recent work by the authors and their colleagues on the development of a second-order 
closure model for high-speed compressible flows is reviewed. This turbulence closure is based 
on the solution of modeled transport equations for the Favre-averaged Reynolds stress tensor 
and the solenoidal part of the turbulent dissipation rate. A new model for the compressible 
dissipation is used along with traditional gradient transport models for the Reynolds heat flux 
and mass flux terms. Consistent with simple asymptotic analyses, the deviatoric part of the 
remaining higher-order correlations in the Reynolds stress transport equation are modeled 
by a variable density extension of the newest incompressible models. The resulting second- 
order closure model is tested in a variety of compressible turbulent flows which include the 
decay of isotropic turbulence, homogeneous shear flow, the supersonic mixing layer, and the 
supersonic flat-plate turbulent boundary layer. Comparisons between the model predictions 
and the results of physical and numerical experiments are quite encouraging. 


x ThiB research was supported by the National Aeronautics and Space Administration under NASA Con- 
tract No. NAS 1-18605 while the authors were in residence at the Institute for Computer Applications in 
Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 
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1. Introduction 


The ability to accurately predict high-speed compressible turbulent flows would have a 
variety of important technological applications in the design of advanced supersonic and hy- 
personic aircraft. Due to the wide range of scales exhibited in these flows, direct numerical 
simulations - with all scales resolved - are impossible for the foreseeable future. Further- 
more, major operational problems with large-eddy simulations in complex wall-bounded 
geometries make their future application to these flows equally questionable. Consequently, 
the method of choice for the calculation of such complex aerodynamic flows continues to 
be based on Reynolds stress modeling and this is likely to remain true for the next several 
decades. Unfortunately, many of the improvements in the Reynolds stress modeling of in- 
compressible turbulent flows that have been developed during the past two decades have 
not, for the most part, found their way into the calculation of the supersonic turbulent flows 
of aerodynamic importance. Design calculations still tend to be based on variable density 
extensions of simple eddy viscosity models such as that due to Baldwin and Lomax [1]. 
However, recent experience with incompressible turbulent flows indicate that significantly 
improved predictions can be obtained from more sophisticated two-equation models and 
second-order closures - particularly in flows with streamline curvature, a system rotation, 
and buoyancy or other body force effects (see Refs. [2 - 4]). While there have been some 
applications of second-order closure models to compressible turbulent flows [5 - 7], these 
studies did not make use of the newer incompressible models and did not explicitly account 
for dilatational effects in a physically consistent manner. This establishes the motivation for 
the present paper which is to provide a much more systematic approach to the development 
of compressible second-order closure models suitable for supersonic turbulent flows. 

The compressible second-order closure model developed in this paper is based on the 
Favre- averaged Reynolds stress transport equation. It will be supplemented with a transport 
equation for the solenoidal part of the dissipation; the compressible dissipation is treated 
using the Sarkar, et al. [8] model which algebraically relates this term to the product of 
the solenoidal dissipation and the square of the turbulence Mach number. Gradient trans- 
port models are used for the Reynolds heat flux, mass flux, and turbulent transport terms 
along traditional lines. However, a new model for the pressure-strain correlation will be 
implemented based on the model developed recently by Speziale, Sarkar, and Gatski [9]. In 
addition, some of the newer developments in near wall turbulence modeling will be made use 
of [10]. An overview of the model’s performance in a variety of test cases will be provided 
along with a brief discussion of needed future directions of research. 

2. The General Model 

Our analysis will be based on the governing equations of motion for an ideal gas which, 
neglecting body forces and bulk viscosity effects, are as follows: 

Continuity 

= 0 ( 1 ) 
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Momentum 
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given that p is the density, u, is the velocity, p is the thermodynamic pressure, R is the ideal 
gas constant, p is the dynamic viscosity, C v is the specific heat at constant volume, T is 
the absolute temperature, cy, is the viscous stress tensor, $ is the viscous dissipation, and 
k is the thermal conductivity. In (1) - (5), the Einstein summation convention applies to 
repeated indices and (•),< denotes a gradient with respect to the spatial coordinate ay. Any 
flow variable T can be decomposed into ensemble mean and fluctuating parts as follows: 

where, for a homogeneous turbulence, the mean J~ can be taken to be a spatial average or, 
for a statistically steady turbulence, it can be taken to be a time average. An alternative 
decomposition based on mass weighted averages can be used wherein 

T = T + F' (7) 

given that T is the Favre average which is defined as 



( 8 ) 


As in the traditional studies of compressible Reynolds stress modeling, both (6) and (7) will 
be used. 

A direct averaging of Eqs. (1) - (3) yields the mean continuity, momentum and energy 
equations which are as follows: 

^ + (pu t ),i=0 (9) 
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where u" = —p'u'Jp is the turbulent mass flux, p'u\ t is the pressure- dilatation correlation, 
and 


are the Favre- averaged Reynolds stress tensor and Reynolds heat flux. Of course, p — pRT 
is the mean pressure. The only assumption made in deriving (9) - (11) is that turbulent 
fluctuations in the specific heat C v can be neglected. 

For high-Reynolds-number turbulent flows, the molecular diffusion terms containing \i 
and k are dominated by the turbulent transport terms except in a thin sublayer near the wall. 
If it is assumed that turbulent fluctuations in the viscosity, thermal conductivity and density 
can be neglected in this region, then the molecular diffusion terms can be approximated as 

follows: 2 

3 (14) 
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Eqs. (14) - (15) are made use of in virtually all of the existing compressible Reynolds stress 
models. The mean viscous dissipation $ can be decomposed as follows: 




(16) 


where e = cr- ii; Jo is the turbulent dissipation rate. Hence, in order to achieve closure, we 

V7 l iJ' r , 

need models for the following turbulence correlations: 


(i) the Favre-averaged Reynolds stress 

(ii) the Favre-averaged Reynolds heat flux Qi 

(iii) the turbulent mass flux u" 

(iv) the turbulent dissipation rate £ 

(v) the pressure-dilatation correlation p'u'.i- 

This is subject to the primary assumption that turbulent fluctuations in the viscosity, thermal 

conductivity, and specific heat can be neglected. 

The Favre-averaged Reynolds stress tensor r \j is a solution of the transport equation 

Q 

+ (P«* T u)j*= -P T ikUj,k ~~P r ik u i,k 
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are the third-order diffusion correlation, the deviatoric part of the pressure-gradient velocity 
correlation, and the dissipation rate tensor. The viscous diffusion term is approximated by 
the relationship 

« ~P( r *j,k + Tikj + r ik%i ) (21) 

which is obtained by an asymptotic near wall analysis neglecting fluctuations in the vis- 
cosity as well as some higher-order correlations, particularly those involving dilatational 
effects. (In a high-Reynolds-number turbulent flow, the viscous diffusion effects would only 
be significant extremely close to the wall where the turbulence Mach number is small and, 
hence, Morkovin’s hypothesis may be invoked.) The pressure-diffusion correlation p f u[ can 
be written in terms of the Reynolds heat flux and mass flux as follows: 


p f u[ = —pRTu" + pRu^T". 

Gradient transport models are used for these quantities in the following form [11] 


( 22 ) 
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where K = is the turbulent kinetic energy and C M = 0.09, a T = 0.7, and a p = 0.5 are 
dimensionless constants. The triple velocity correlation is modeled by using a traditional 
gradient transport hypothesis as follows [4] 


2 2 

= ~1 C >— fa;.* + T **.; + *-;*..•) (25) 

where C, = 0.11. Eq. (25) is the isotropized version of the gradient transport model 
introduced by Launder, Reece, and Rodi [12]. The advantage of this formulation lies in the 
fact that each turbulent transport term is coupled with a viscous term of the same form — 
a feature that arises from the isotropization of the model for the triple velocity correlation 
and the use of the true dissipation rather than the pseudo-dissipation in the formulation of 
the viscous terms in (21). 

Now, we have to discuss the modeling of the remaining correlations needed for closure 
which consist of the deviatori c par t of the pressure gradient- velocity correlation II;j, the 
pressure-dilatation correlation and the dissipation rate tensor Since our model is 

developed for high turbulence Reynolds numbers, the Kolmogorov assumption of isotropy of 
dissipation will be invoked sufficiently far from solid boundaries, i.e., 

2_ c 

P^j = g pzoij (26) 
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where ~pe = is the scalar dissipation rate of the turbulence which appeared in Eq. 

(16). Following the work of Sarkar et al. [8], the dissipation rate will be decomposed into 
solenoidal and compressible parts as follows: 

£ = c a £ c (27) 

where, for a homogeneous turbulence, = pu)[u)[ and ~p£ c = |/i(u^) 2 given that is the 
fluctuating vorticity. A comparable decomposition was recently introduced by Zeman [13] 
in an analysis of eddy shocklets. Based on a simple asymptotic analysis, a model for the 
compressible dissipation was derived by Sarkar et al. [8] wherein £ c = cl\M 2 s m given that 
Mt = (2K j^RT) 1 / 2 is the turbulence Mach number, 7 = (7 P /Z7 V is the ratio of specific heats, 
and ai is a dimensionless constant. A calibration based on direct numerical simulations of 
compressible isotropic turbulence yielded the final model [8] 

e = (l + a a M 2 )£ 8 (28) 


with a.\ = L The asymptotic analysis that led to the development of (28) also indicated that 
the compressible dissipation is substantially larger than the pressure dilation correlation in 
compressible isotropic turbulence [8]. Consequently, at this stage of the modeling we neglect 
the pressure dilatation correlation, i.e., we set 


p' u i,i 


0. 


(29) 


Work is currently underway to develop more general models for the pressure-dilatation cor- 
relation suitable for high-speed flows with strong mean velocity gradients. 

The exact transport equation for the solenoidal dissipation is of the general form 


d" (/?&»£, )„* — ^Pe t H" T^e 9 d“ 


(30) 


given that fluctuations in the viscosity can be neglected. Here, V e , is the production of 
dissipation due to deviatoric strains, $ ej is the destruction of dissipation, and T> Cl is the 
turbulent transport of dissipation (the precise form of the expressions for V e ,,^ e , and T> Ct are 
quite complicated and are omitted for simplicity). Consistent with the asymptotic analysis 
for the solenoidal dissipation, the higher-order correlations in (30) will be approximated by 
a variable density extension of the commonly used incompressible models. This yields a 
modeled transport equation for e, of the form 


d(/>£«) 

dt 


14 £ 2 

~ “ @c2P 


(31) 


H“(CeP 

£a 

where C e \ = 1.44, C e 2 = 1.83, and C c = 0.15 are dimensionless constants. In the incom- 
pressible limit, the standard form of the modeled dissipation rate equation is recovered [12]. 
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It should be noted that, unlike in most previously proposed compressible models, mean di- 
latational effects are accounted for exactly by the second term on the r.h.s. of (31) - a feature 
that allows for the description of compressed isotropic turbulence as we will see later. 

Finally, in order to close these equations, a model for the deviatoric part of the pressure 
gradient-velocity correlation Ily is needed. We will assume that, to the first-order, II, ^ can 
be approximated by a variable density extension of its incompressible form (i.e., the leading 
order effects of compressibility are through changes in the density). We will use a variable 
density extension of the SSG model which was recently derived by Speziale, Sarkar, and 
Gatski [9] for incompressible turbulent flows. This model takes the form: 

n*j = p[(Ci£ + C*7*)bij — C , 2e(fcifcl>jfcj — — 


—CzK(hikSjk + bjkSik 


g ^mn *Sfnn bij ) 


-C 4 K(b ik W jk + b jh Wi k ) 


(32) 


where 


-f 1 - C’// I/s )(S., - I&A,)] 

(t» - 

,J 2 K 

$ij — 2 (^»«i %.*)j = **•*"■* ) 

V = -TijUij, II = hjbij 


(33) 

(34) 

(35) 


and Ci = 3.4., C[ = 1.80, C 2 = 4.2, C 3 = 1.25, C 4 = 0.40, C* - 1.62 are dimensionless 
constants. The SSG model was tested in a variety of incompressible, homogeneous turbulent 
flows where it was shown to yield improved predictions over older models such as the Launder, 
Reece, and Rodi model (see Speziale, Sarkar, and Gatski [9] for more details). 

The complete form of our proposed compressible second-order closure model has now 
been provided for high turbulence Reynolds number flows. However, some discussion is 
needed concerning the integration of this model to a solid boundary. For attached boundary 
layers, wall functions can be used to bridge the outer and inner flows. If this is done, some of 
the newer wall functions that have been developed for supersonic turbulent boundary layers 
should be made use of (see He, Kazakia, and Walker [14] for an interesting discussion). On 
the other hand, with minor modifications, this second-order closure model can be integrated 
directly to a solid boundary - the preferred approach when there is separation or other 
complex alterations of the wall boundary conditions. In order to integrate the second-order 
closure model presented in this study to a solid boundary, two general alterations to the 
model are needed: 


(a) the destruction of dissipation term in the modeled dissipation rate equation (31) must 
be damped to zero as the wall is approached. More precisely, C e2 must be replaced by C t2 f 2 
where f 2 — » 1 for Re t 1 and f 2 — ► 0 as y — > 0 (here, y is the coordinate normal to the wall 
and f 2 is 0(y 2 ) near the wall). 
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(b) A near wall correction to the deviatoric parts of the dissipation rate tensor and 
pressure gradient-velocity correlation must be provided. 

Variable density extensions of existing incompressible models can be used as a first approx- 
imation (the reader is referred to Speziale, Abid, and Anderson [10], Launder and Shima 
[15], and Gatski et al. [16] for a more detailed discussion of near wall modeling). However, 
these near wall models involve wall damping functions that tend to be quite ad hoc. Con- 
sequently, a new approach is currently under development wherein only the momentum and 
turbulent kinetic energy equations are integrated to the wall; the individual components of 
the Reynolds stress tensor, as well as the dissipation rate equation, are only integrated to 
the edge of the viscous sublayer eliminating most of the need for wall damping. 

Finally, it should be noted that for many practical applications to supersonic flows, it is 
more convenient to integrate the total energy equation for E = C V T + \uiU{ rather than the 
thermal energy equation (11). This equation takes the form 


Q 

oftpE) + (7*W5)>i = + W ‘J U j - PUi 

- 9i),i +( CT ij U 'j - p' u i - P E " U ")n 
where the energy flux is given by 

+ TijUj + 


(36) 


(37) 


and the other correlations in (36) are modeled as previously described. 


3. Applications of the Model 

The first application to be considered is the decay of compressible isotropic turbulence 
with no mean dilatation (i.e., all components of Uij are zero). In Figures 1(a) - (b), the 
results predicted by the model for the decay of the turbulent kinetic energy are compared 
with those obtained from direct numerical simulations (DNS) for initial turbulence Mach 
numbers Mt , o of 0.1, 0.3, and 0.4 (these results are taken from Sarkar et al. [8]). It is clear 
that the model does an excellent job in reproducing the trends of the DNS which indicate 
a moderate increase in the decay rate of the turbulent kinetic energy with increasing turbu- 
lence Mach number. (The results should not be compared from a quantitative standpoint 
since the model is for high Reynolds number turbulence and the DNS is for R\ = 15). This 
Mach number dependence of the turbulence decay rate is due to the compressible dissipa- 
tion and cannot be described by the commonly used compressible turbulence models which 
neglect such dilatational effects. In Figure 2, the time evolution of the pressure-dilatation 
correlation obtained from the DNS of isotropic turbulence [8] is shown for an initial turbu- 
lence Mach number of M tfi = 0.5. It is clear that after the early transient effects due to the 
initial conditions die out, the pressure-dilatation correlation becomes extremely small. This 
validates the assumption made of neglecting the pressure-dilatation to the lowest (isotropic) 
order. 
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The second problem that we will consider is the rapid compression or expansion of an 
isotropic turbulence. In this problem, a decaying isotropic turbulence is subjected to the 
mean dilatation 




-T 
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0 

0 

\ 

II 

• •• 


0 

ir 
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0 




0 

0 

-T 

3 

/ 


(38) 


where T is a constant (T > 0 corresponds to an expansion whereas T < 0 represents a 
compression). We will consider rapid compressions and expansions for which \T\Ko/eo 1. 
The general model reduces to the following simple differential equations for this problem: 


K = - fa - (1 + ai M})e. 


(39) 


£. = 



(40) 


For |r|if 0 /eo > 1) the short-time solution to Eqs. (39) - (40) can be approximated as 


K = Koexp(-^Tt) (41) 

4 

= £ a , 0 exp(--ri) (42) 


since the first term on the r.h.s. of (39) and (40) dominate. Hence, the integral length scale 
predicted by the model is as follows: 


A = A 0 exp 



(43) 


since A oc K 3 ^ 2 /e a . 

The results given in (41) and (43) are identical to those obtained by Reynolds [17] based 
on Rapid Distortion Theory (RDT). As correctly pointed out by Reynolds [17], many of the 
commonly used compressible turbulence models, unlike our model, erroneously predict that 
the integral length scale decreases with an expansion and increases with a compression. This 
problem arises from the neglect of the exact mean dilatational term — |pe 4 u, i that appears 
on the r.h.s. of (30). 

Now, we will consider the problem of compressible, homogeneous shear flow. The model 
predictions will be compared with the direct numerical simulations conducted recently by 
Sarkar, Erlebacher, and Hussaini [18]. In this problem, an initially decaying compressible 
isotropic turbulence is subjected to a uniform shear S with the corresponding mean velocity 
gradients 


u i,j = 


0 S 0 \ 
0 0 0 
0 0 0 / 


(44) 


In Figure 3, the time evolution of the turbulent kinetic energy predicted by our second- 
order closure model is compared with the predictions of the Launder, Reece, and Rodi 
(LRR) model and the results of direct numerical simulations [18] for the initial conditions 
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SK 0 /e = 7.18, R\ 0 = 15, and M to = 0.2. It is clear from this figure that the new model 
does a better job in reproducing the trends of the simulations. (Some differences would be 
expected due to the fact that the model is for high turbulence Reynolds numbers whereas 
the simulations are for low R\). As discussed earlier, the compressible dissipation model 
(28) is based on the scaling 

— ^ — = constant. (45) 

Mfe, 

In Figure 4, the DNS results of Sarkar et al. [18] for the time evolution of eJMfe, is shown 
for a variety of initial conditions. After the early time transient, e c /M?e, settles down to a 
constant value in the range of 0.5 to 0.6 thus validating the scaling assumed in (45). (We still 
leave a.\ = 1 to account for the fact that we have neglected the pressure-dilatation correlation 
- a term which in shear flow has a comparable effect to the compressible dissipation.) 

The next problem that we will consider is the plane supersonic mixing layer. In this 
problem, a supersonic stream mixes with a subsonic stream downstream of a splitter plate 
(see Figure 5). The calculations that we will present are an extension of those conducted 
by Sarkar and Balakrishnan [19] for a high-speed stream velocity U x = 2500 m/sec and a 
low-speed stream velocity U 2 = 800 m/sec. Both streams have the same thermodynamic 
quantities which are prescribed as follows: T\ — T 2 = 800°#, pi = Pa = 1 atm, and 
Pi — Pi = 0-44 kg/m 3 . A convective Mach number M c and spreading rate C$ are defined by 


Pi - u, 

M c = 

(Xi + &2 

(46) 

d8 „ fU\ — U 2 \ 
di ~ S \lh + Uj 

(47) 


where 8 is the thickness of the shear layer and a is the local speed of sound. In Figure 
6, the normalized spreading rate C&/(Cs) o obtained from our second-order closure model is 
compared with the experimental “Langley curve” for a variety of convective Mach numbers 
M c . The comparison between the model predictions and experiments is quite good. In 
contrast to these results, a variable density extension of either the Launder, Reece, and Rodi 
model or the SSG model were not able to predict the dramatic drop in the spreading rate 
that occurs for convective Mach numbers M c > 1 as shown in Figure 6. Consequently, it is 
clear that dilatational effects must be directly accounted for if supersonic shear layers are 
to be properly described. In Figures 7 - 8, the normalized streamwise turbulence intensity 
and turbulent shear stress predicted by the second-order closure model are compared with 
experimental data (here, AC7 =Ui~ U 2 ). The predictions of the Launder, Reece, and Rodi 
model are also displayed. It is clear from these results that our model predictions are well 
within the range of the experimental data unlike those of the LRR model. 

Finally, we will present some preliminary results for the supersonic flat plate turbulent 
boundary with an adiabatic wall and with wall cooling [20]. These calculations were done by 
a direct integration of this second-order closure model to the wall using damping functions 
as discussed at the end of Section 2. In Figure 9, the normalized skin friction predicted 
by the second-order closure model is compared with the Van Driest curve and experimental 
data for a range of external Mach numbers 0 < M e < 10 (the Reynolds number based on 
the momentum thickness Re fi = 10 4 ). Clearly, the model predictions are in good agreement 
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with the experimental data. In Figure 10, the normalized skin friction predicted by this 
second-order closure model for the supersonic fiat plate boundary layer with wall cooling is 
shown (M e - 4, Re e = 10 4 , T w = wall temperature, and T aw = adiabatic wall temperature). 
Again, the model predictions are in reasonably good agreement with the Van Driest curve. 
More detailed calculations of supersonic turbulent boundary layers are currently underway 
(Speziale and Abid [20]). 

4. Concluding Remarks 

A new and reasonably comprehensive second-order closure model has been developed for 
high-speed compressible turbulent flows. The features that distinguish this new model from 
Previously proposed second-order closures for compressible turbulence lies in the explicit 
treatment of dilatational effects on the turbulence dissipation rate and the use of recently 
improved models for the pressure gradient-velocity and viscous diffusion correlations. Pre- 
liminary tests of the model in a variety of homogeneous and inhomogeneous supersonic 
turbulent flows were quite encouraging. 

A substantial amount of research is still needed to refine the model. The major issues 
that need to be further explored are as follows: 

(1) the development of an anisotropic correction to the model for the pressure-dilatation 
correlation valid for significant mean strain rates, 

(2) the development of modeled transport equations for the Reynolds heat flux and mass 
flux terms that transcend the simple gradient transport models now being used, and 

(3) the development of better near wall models that incorporate more compressible turbu- 
lence physics and are numerically robust. 


In addition, the ability of this second-order closure model to predict the turbulence 
statistics in compressible flows with strong shocks remains to be demonstrated. Furthermore, 
if this model is to be extended to the hypersonic flow regime, real gas effects have to be 
incorporated into the model. All of these issues are currently under investigation in our 
ongoing research effort on compressible turbulence modeling. 
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Figure 2. Time evolution of the pressure-dilatation correlation in compressible isotropic 
turbulence obtained from Direct Numerical Simulations [8] (R\o = 15, M to = 0.5). 
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Figure 3. Time evolution of the turbulent kinetic energy in homogeneous shear flow (R\ 0 = 

15, Mt 0 = 0.2, SKo/eo = 7.18): Direct Numerical Simulations [18]; - • - LRR Model; 

• • • Current Model. 
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Figure 4. Time evolution of eJMfe, in homogeneous shear flow obtained from Direct Numer- 
ical Simulations [18] for a variety of initial conditions (the initial turbulence Mach number 
M to varies from 0.2 to 0.5). 
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Figure 5. Schematic diagram of the supersonic mixing layer. 
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Figure 6. Normalized spreading rate as a function of the convective Mach number M c \ 
comparison of the model predictions with experimental data for the supersonic mixing layer. 
Current Model; LRR Model; a Experimental “Langley Curve” (Refs. [16,19]). 
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maximum streamwise turbulence intensity as a function of the con- 
Current Model; LRR Model; a Experimental Data (Refs. 
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Figure 8. Normalised maximum turbulent shear stress as a function of the convective Mach 
m er. urrent Model; - - LRR Model; a Experimental Data (Refs. [16,19]). 
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Figure 9. Normalized skin friction as a function of the external Mach number M e for the 
compressible flat plate turbulent boundary layer with an adiabatic wall (Re$ = 10 4 ): o Model 
Predictions; Van Driest Curve. 
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Effect of Surface Cooling 
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Figure 10. Normalized skin friction as a function of wall temperature T w /T aw : the effect of 
surface cooling on the compressible flat plate turbulent boundary layer for an external Mach 
number M e — 4 (Reg = 10 4 ). o Model Predictions; Van Driest Curve. 
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